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^  ABSTRACT  j 

A  collocation  method  is  given  for  steady-state  simulation  of 
multiple  reactions  in  porous  catalysts.  A  realistic  multicomponent 
diffusion  model  is  used,  which  includes  an  allowance  for  pore  size 
distribution.  Hyperbolic  basis  functions  are  introduced  to  repre¬ 
sent  the  intraparticle  profiles;  compact  solutions  are  thus  obtained 
both  in  the  presence  and  absence  of  fast  reactions.  Calculations 
for  a  six-component  catalytic  reforming  system  show  that  the  catalyst 
performance  is  strongly  affected  by  intraparticle  diffusion. 
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SIGNIFICANCE  AM)  EXPLANATION 


Theoretical  models  of  gas  diffusion  and  flow  in  porous  solids 
are  well  developed,  and  are  beginning  to  be  applied  to  simple 
catalytic  systems.  Calculations  of  this  kind  permit  a  new  level 
of  mderstanding  of  catalysis,  which  should  lead  to  more  efficient 
chemical  processes. 

Detailed  simulations  of  catalyst  particles  should  be  especially 
useful  in  studies  of  multi-reaction  processes,  for  which  the  catalyst 
selectivity  may  be  sensitive  to  intraparticle  diffusion.  In  this 
paper  we  suimarize  the  relevant  transport  equations  and  give  a  new 
collocation  method  for  solving  this  kind  of  problem. 

The  simulation  of  a  catalyst  particle  needs  to  be  done  effi¬ 
ciently  if  it  is  to  be  included  in  a  reactor  computation.  Colloca¬ 
tion  methods  based  on  global  polynomials  become  inefficient  in  the 
presence  of  fast  reactions  because  of  the  steep  reaction  fronts 
which  then  occur.  Inproved  basis  functions  are  introduced  here 
from  a  linearized  version  '  *  the  transport  equations,  this  pro¬ 
viding  compact  solutions  .  ,/ide  range  of  reaction  rates. 


The  responsibility  for  the  wording  and  views  expressed  in  this 
descriptive  surma ry  lies  with  MFC,  and  not  with  the  authors  of  this 
report. 
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COLLOCATION  ANALYSIS  OP  MULTICOMPONENT  DIFFUSION 
AND  REACTIONS  IN  POROUS  CATALYSTS 


Jan  P.  SArensen  and  Warren  E.  Stewart* 


INTRODUCTION 

Theoretical  models  of  gas  diffusion  and  flow  in  porous  solids  are 
well  developed  [23,  24,  10,  11,  17],  and  are  beginning  to  be  applied 
to  simple  catalytic  systems  [2,  9,  15,  18,  19,  20,  38].  Calculations 
of  this  kind  permit  a  new  level  of  understanding  of  catalysis,  which 
should  lead  to  more  efficient  chemical  processes. 

Detailed  simulations  of  catalyst  particles  should  be  especially 
useful  in  studies  of  multi-reaction  processes,  for  which  the  catalyst 
selectivity  may  be  sensitive  to  intraparticle  diffusion.  In  this  paper 
we  summarize  the  relevant  transport  equations  and  give  a  new 
collocation  method  for  solving  this  kind  of  problem. 
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The  simulation  of  a  catalyst  particle  needs  to  be  done 
efficiently  if  it  is  to  be  included  in  a  reactor  computation. 
Collocation  methods  based  on  global  polynomials  become  inefficient  in 
the  presence  of  fast  reactions  [35*  26,  7,  36]  because  of  the  steep 
reaction  fronts  which  then  occur.  Improved  basis  functions  are 
introduced  here  from  a  linearized  version  of  the  transport  equations, 
thus  providing  compact  solutions  over  a  wide  range  of  reaction  rates. 

TRANSPORT  EQUATIONS 

Consider  the  steady  diffusion  and  reaction  of  a  multicomponent 
gaseous  mixture  in  a  porous  catalyst  particle.  We  use  the  following 
model  [24,  10,  11]  to  describe  the  intraparticle  molar  fluxes: 

■v 

»<">  .  -  Yjnj W  [Pair's*  -  (B0/m)  -  Da£o  (1  ) 

lc=1 

The  first  right-hand  term  arises  from  gaseous  diffusion,  the  second 
from  viscous  flow,  and  the  third  from  surface  diffusion.  The  gaseous 
diffusion  expression  includes  the  leading  thermal  transpiration  term 
n^airtT/az  from  the  dusty-gas  model  of  Mason  et  al.  [23].  The  elements 
of  F(r  )  are  obtained  by  setting  r  =  r  in  the  expressions 

K  K 

nc 

Pu(r)  •  1/VD  -  IV-t  (2a) 

h=1 

h?<i  i  =  1  ,...n 

c 

j  =  1 , . .  «n 

C 

Pij(r)  “  “  ci/cBij  1  *  J  (2b) 

in  which  n^  is  the  mmber  of  gaseous  species. 
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Equation  (1 )  includes  the  Knudsen  diffusion  equation  and  the 
Stefan-Maxwell  equation  [6,  16]  as  asymptotes  for  low  and  hi$i 
pressures,  respectively.  This  model  has  also  been  tested  thoroughly 
over  the  intermediate  region  [10,  11 ]  commonly  encountered  in 
industrial  gas-solid  processes.  The  summation  on  k  in  eqn  (1 )  is  a 
quadrature  of  the  flux  expressions  of  Mason  et  al.  [24]  over  the  pore 
size  distribution.  A  two-point  sum  (n^  =2)  fits  the  available  data 
well  for  a  variety  of  catalysts  [10,  11  ];  one  point  suffices  for 
narrow  pore  size  distributions. 

Equation  (1 )  was  developed  for  non-reacting  systems,  under  the 
assumption  that  the  concentrations  and  temperatures  in  the  pores 
conform  to  smooth  fields  defined  throughout  the  porous  medium.  For 
simplicity  we  use  the  same  model  in  the  presence  of  chemical 
reactions.  This  approach  should  be  satisfactory  if  the  pores  are 
cross-linked  frequently,  so  that  the  concentration  changes  along 
individual  pore  segments  are  small.  This  point  is  discussed  in  greater 
detail  by  Jackson  [17];  see  also  Mingle  and  Smith  [25]. 

The  energy  flux  within  the  particle  is  modelled  as 

nc 

*  r  ?i  ®i  -  v« 21  °> 

i=1 

that  is,  as  a  sum  of  convective  and  conductive  terms,  with  the  Dufour 
effect  [6]  neglected. 

In  the  following  calculations,  we  use  the  interstitial 
concentrations  c^  and  the  temperature  T  as  state  variables.  Equations 
(1 )  and  (3)  are  readily  rewritten  in  these  variables  hy  insertion  of 
the  ideal  gas  law.  The  result  can  be  expressed  in  the  matrix  form 


The  flux  array  N  has  the  elements  (N  ,  ...  N  ,  N  j ,  and  the  state 
~  ~1  ~nc  ~E 

vector  £  has  elements  |c. ,  ...  c  ,  T} .  The  elements  of  D  are 

1  nc  ~ 

summarized  in  Appendix  A. 

At  steady  state,  the  fluxes  in  the  particle  satisfy  the  mass  and 
energy  balances 


in  the  smooth-field  approximation.  Here  is  the  local  production 

rate  of  species  i  per  unit  volume  of  the  heterogeneous  median;  its 

evaluation  for  fairly  general  kinetics  is  discussed  in  [33]* 

2  p 

In  this  work,  we  consider  symmetric  states  {c^Cz  ),  T(z  )}  in  a 

catalyst  slab  or  sphere,  with  boundary  conditions 

c  1  =  c,.;  dc Vdz  I  =0  i*1,...n  (7) 

4Ul  10  1  I  z=0  c 

T  I  =  T  ;  dT/dz  =  0  (8) 

j  z=1  z=0 

The  dimensionless  coordinate  z  is  measured  from  the  particle  center, 
as  a  fraction  of  the  particle  half-thickness  L.  More  general  boundary 
conditions  are  readily  accomodated. 

A  collocation  procedure  for  solving  eqns  (4)-(8)  is  described  in 
the  next  four  sections.  After  the  description,  three  numerical 
examples  are  given. 
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BASIS  FUNCTIONS 


The  collocation  procedure  requires  a  set  of  basis  functions  to 
describe  the  concentration  and  temperature  profiles.  Global 
polynomials  are  often  used  for  this  purpose  [35,  12,  36],  but  many 
terms  are  then  required  to  get  acceptable  accuracy  for  fast  reactions. 
Piecewise  polynomials  (splines)  are  more  flexible  [26,  7,  36],  but  the 
choice  of  breakpoints  for  these  functions  needs  further  study, 
especially  for  multicomponent  problems.  The  approach  taken  here  is  to 
develop  natural  basis  functions  from  a  linearized  form  of  the  given 


problem.  A  similar  approach  has  been  applied  successfully  to  systems 
of  first-order  differential  equations  [13]. 

Insertion  of  eqn  (4)  into  (5)  and  (6),  and  linearization  around  a 
reference  state  6^,  gives  the  matrix  differential  equation 


D(S  ) 

~  ~R 


y2s 


=  R(£  ) 
~  ~R 


+  R 


(5  ) 

~R 


[  ^  - 


5  ] 

~R 


(9) 


Here  D(^  )  is  the  transport  coefficient  matrix  calculated  at  S  from 
~  ~R  ~R 

the  relations  in  Appendix  A,  R(£  )  is  the  vector  of  production  rates 

~  ~R 

at  f  ,  and  R*  ( £  )  is  the  matrix  9R/3£  evaluated  at  £  from  the  given 
n  ~  ~R  ~  ~  ~R 

kinetic  model.  The  matrix  D(£  )  is  non-singular,  since  eqns  (l)  and 

~  ~R 


(3)  are  linearly  independent  [31  ]. 

The  solution  vector  |  of  eqns  (7)-(9)  is  a  linear  combination  of 
the  following  functions: 


1 ,  cosh(z\/""x^) 

1,  sinh(z\/"“x^}/z 


(Slab)  k  =  1 ,  ...  r 

(10) 

(Sphere)  k  =  1 ,  ...  r 


Here  r  is  the  rank  of  the  matrix  [-  D(£  )-1  L2  R'U  )  ],  and 

~  ~R  ~  ~R 

y  •••  Xf  are  its  nonzero  eigenvalues.  If  a  p-fold  eigenvalue  yields 
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fewer  12ian  p  eigenvectors,  alternative  functions  given  in  eqn  (11a) 
or  (11b)  will  appear.  Since  in  practice  the  rank  may  be  uncertain,  we 
use  its  upper  bound  m:  the  maximum  number  of  independent  production 
rates  permitted  by  the  stoichiometry  and  local  constraints  of  the 
reaction  system  [3,  16,  17,  31,  33]*  In  the  stoichiometric  analysis 
of  [33],  m  is  the  number  of  non-equilibrium  reactions  that  yield  pivot 
coefficients  for  mobile  species.  The  m  largest  values  of  Real(\/  \  ) 

A 

are  then  used  to  obtain  m  functions  of  the  form  in  eqn  (10). 

The  quantities  V~X^  axe  generalized  Thiele  moduli  [4]  for  the 
corresponding  eigenfunctions  in  eqn  (10).  Polynomial  approximations  of 
the  linearized  solution  became  difficult  if  any  of  these  moduli  are 
large;  in  such  cases  we  will  say  that  the  differential  equations  are 
"stiff". 

Equation  (10)  is  written  for  non-zero  eigenvalues  X  ,  ...  X  , 

^  m 

i.e.,  for  kinetics  of  full  rank  m.  In  this  case,  the  spatial  function 

"1"  corresponds  to  an  equilibrium  solution  £  of  eqn  (9)-  If  any  of 

,  ...  in  eqn  (10)  are  zero,  then  the  corresponding  functions  are 

replaced  by  z^,  obtained  by  expansion  of  the  solutions  in  powers  of 

2 

\f\.  This  procedure  provides  a  particular  solution  proportional  to  z 
whenever  eqn  (9)  lacks  an  equilibrium  solution;  in  such  cases,  the 
function  "1 "  is  merely  a  solution  of  the  related  homogeneous 
equation,  ^7  £  =  0.  Kinetics  of  rank  less  than  m  can  arise  from 
zero-order  rate  expressions  (which  we  prefer  to  avoid) ,  or  from  other 
causes  such  as  absence  of  various  concentrations  from  the  reaction 
rate  expressions. 

For  non-linear  problems,  we  extend  the  function  set  of  eqn  (10) 

by  differentiation  with  respect  to  \/\  at  each  characteristic  value 

X  •  The  following  differential  forms  are  convenient: 

1c 
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1 - cosh(zWT) 

,  avr/ 


(1-z)s  (1+z)S 

=  -  exp(zVxT)  + -  exp(-zWTT ) 

Xk  2  2 


3  —  1  f  2|  •  •  • 


(Slab)  (11a) 


sinh(z\ZT) 

z 


(1-z)3 

-  exp(z\/T7) 

2z  K 


(1«)S  __ 

-  exp(-zyx7) 

2z  K 


s  =  1,  2,  ...  (Sphere)  (11b) 

These  functions  are  also  used  to  deal  with  clustered  eigenvalues  as 

described  below.  The  solutions  E  (z)  and  N  (z)  are  then  approximated 

i  iz 

as  combinations  of  basis  functions  selected  from  eqns  (10)  and  (11), 


n 

a  $  (z)  i  ■  1,  ...  nc+1  (12) 

j=0 
n 

"iz  -  E  fu  fc  1  =  1.  •••  I>c+1  (13) 

j=1 

with  adjustable  coefficients  a^  and  f ,  and  with  $Q(z)  =  1.  The 
index  n  is  chosen  by  the  user.  Here  and  below,  we  mark  the  approximate 
solutions  with  a  tilde  (~).  Equation  (13)  is  included  to  facilitate 
the  treatment  of  eqns  (5)  and  (6)  for  systems  with  variable  D. 

For  sufficiently  large  a,  eqn  (12)  can  approximate  any  continuous 

p 

symmetric  function  C  (z  )  to  arbitrary  accuracy  over  the  interval 
[0,1],  even  if  all  the  values  \J  are  replaced  by  a  single  constant 
a.  To  prove  this,  we  note  that  eqn  (12)  then  reduces  to  a  repre¬ 
sentation  of  -  aiQ) [Slab]  or  -  aiQ) [Sphere]  by  an 
expansion  P^O  -z)exp(az)  +  Pn(l+z)exp(-oe).  Division  by  exp(az)  leads, 
for  either  geometry,  to  a  representation  of  a  continuous  function  by  a 
polynomial  P  (1-z);  the  proof  then  follows  directly  from  the 


-8- 


Weieratrass  approximation  theorem  [8].  Equation  (13)  has  similar 
approximating  power  for  the  fluxes  N.^.  Consequently,  it  ib 


permissible  to  modify  the  X's  for  simplicity  if  a  sufficient  number 
of  basis  functions  is  used. 

The  following  ordering  of  the  functions  for  selection  has  given 

good  results,  let  the  desired  number  of  collocation  points  be  n;  then 

m-1  functions  4>^(z)  are  required.  We  start  with  $q(z)  =  1  and  =  0; 

then  we  choose  the  constants  ,  ...  as  the  values  of  Real(\/  X^)  in 

descending  order.  These  values  are  necessarily  non-negative.  If  n  <  m, 

we  drop  the  values  after  ot^,  thus  omitting  those  terms  of  the  linearized 

solution.  On  the  other  hand,  if  n  >  m,  we  insert  a. ,  •••  an  according 

to  the  recursion  formula  a.  =  a.  .  Then  we  rearrange  and  relabel  the 

J  J-® 

resulting  list  to  form  an  ascending  sequence  {a  I  =  {art,  ...  a  }. 

j  0  n 

Close  groups  of  unequal  a  values  are  then  compacted  as  follows, 

to  strengthen  the  linear  independence  of  the  basis  functions.  If  t^y 

a  <  2.0j  for  a  slab,  or  a  <  2.5j  for  a  sphere,  with  o  <  o  ,  then 
J  “  J  ""  '  J 

a  and  all  equal  or  smaller  o's  are  replaced  by  zeros,  thus  replacing 
J 

those  basis  functions  by  polynomials.  Subsequent  sequences  {  a^,  ... 

a  },  with  k  the  largest  integer  such  that  a.  /a  <  (1 .4)k_1 ,  are 
h+k  h+k  h  — 

compacted  by  replacing  each  member  with  (<*h  . . .  ah+k)  1  .  These 

procedures  were  developed  from  numerical  tests  to  provide  a 
well-conditioned  Gram-Schmidt  solution  of  eqn  (19). 

The  basis  functions  are  then  obtained  by  rewriting  eqns  (10)  and 
(11)  as  follows,  with  multipliers  exp(-a  )  to  ensure  computable  values: 

J 
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For  a  =  0  in  a  slab  or  sphere, 
3 

.  =  z  J 

J 


Sj 


(14a) 


For  a  >  0  in  a  slab, 

3 

S  •  s 

(1-z)  ^  (1+z''  ^ 

<j>  _  - - -  exp[a  . ( z— 1 )]  + - - -  exp[-a.(z+1 )]  (14b) 

J  2  0  ^  0 


For  a  >  0  in  a  sphere, 

j 

(1-z)  (1+z) 

<t>  -  -  exp[  a  (z-1 )] - exp[-a  (z+l)J  (14c) 

3  2z  j  2z  j 


Each  integer  s  is  taken  here  here  as  the  number  of  prior  occurrences 
3 

of  the  associated  value  a  in  the  a-list. 

3 


The  selection  order  used  here  is  different  from  that  used  by 
Guertin  et  al.,  who  gave  last  priority  to  rapidly  decaying  functions  in 
their  basis  selection  for  initial  value  problems.  Such  a  rule  is  not 
relevant  for  the  boundary-value  problems  considered  here,  in  which  the 
fluxes  depend  on  the  gradients  and  thus  the  steep  basis  functions 
may  be  important.  By  giving  priority  to  the  large  a's  (after  a^) ,  we 
bracket  any  unused  values  and  obtain  good  accuracy  with  fewer 
basis  functions. 


COLLOCATION  POINTS 

The  collocation  procedure  consists  of  adjusting  eqns  (12)  and 

(13)  to  satisfy  eqns  (7)  and  (8)  exactly,  and  to  satisfy  eqns  (4)-(6) 

at  a  set  of  interior  points  z  ,  . . .  z  .  We  choose  these  points  by 

1  n 

analyzing  the  residuals  of  eqns  (3)  and  (6), 

7.1  =  1-  [  zatf.  (£)]  -  L  R  (£)  i  =  1,  ...  n  +1 

1  “  za  dz  12  1  ~  c 

(15) 


in  the  series  form 
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00 


?l£  -  «n<z>  Z 


i  =  1 ,  . . .  n  +1 
c 


(16) 


k=0 


Here  the  d  are  unknown  coefficients,  and  Q^Cz)  is  a  function  whose 


zeros  z  ,  ...  z  will  he  the  interior  collocation  points.  We  choose 
1  n 


n-1 


Q  (z)  =  <t>  (z)  +  7  h.  <f>  (z) 
n  n  L  ]  3 

j=0 


(17) 


in  which  the  h  are  coefficients  to  be  determined, 

j 


We  minimize  the  magnitudes  of  the  mean  residuals 

'i  r\ 


T  r  za  dz  =  H.  (1 )  -  L 
•'i  »  iz 


R^(|;)  za  dz 


1 


1  oo 


(18) 


Y.  dU£  Vz)  *kU)  iz 


i  =  1 ,  . . .  n  +1 

c 


0  k=0 


with  respect  to  the  coefficients  b^  by  imposing  the  following 


orthogonality  conditions  on  Qr(z) : 


/: 


Q  (z)  (j>,  (z)  za  dz  =  0 
n  k 


k  =  Oj  • • •  ihI 


(19) 


These  conditions  cause  Q  (z)  to  have  n  distinct  zeros  on  the  interval 

n 


(0,1 ),  as  shown  in  Appendix  B.  Collocation  at  the  zeros  of  Q^(z)  thus 
is  feasible,  and  eliminates  the  first  n  terms  of  the  mean  residuals 


as  expressed  in  eqn  (18).  The  coefficients  bj  in  eqn  (17)  are 


determined  by  eqn  (19);  the  values  of  the  d  are  not  required. 

XA 


Equations  (19)  are  also  obtainable  by  minimizing  the  mean  square 


of  Q  (z)  with  respect  to  b.,  ...  b  . .  The  leading  term  of  each 
n  u  n-i 
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residual  in  eqn  (16)  is  thus  minimized  in  a  least-squares  sense 

over  the  interior  of  the  catalyst  particle. 

In  systems  with  constant  D,  eqns  (13)  are  exactly  consistent  with 

eqns  (12)  and  (4).  In  systems  with  variable  D,  however,  flux  residuals 

®i  -  5^(|)  d|/dz  arise.  The  resulting  truncation  error  in  eqn  (20)  is 

normally  small,  and  could  be  reduced  by  using  a  second  set  of  points 

for  evaluation  of  ft  .  This  question  will  be  considered  at  another 
iz 

time. 


DISCRETE-ORDINATE  EQUATIONS 


It  is  convenient  to  compute  the  solution  values  directly  at  the 
collocation  points,  rather  than  solve  for  the  coefficients  in  eqns 
(12)  and  (13).  This  is  done  by  expressing  the  gradient  and  divergence 
operations  in  the  forms 


(20) 


(21) 


with  a  =  0  for  slabs  and  a  =  2  for  spheres.  Equation  (21 )  has  n  terms, 
as  does  eqn  (13);  for  simplicity  the  needed  n  values  of  are  taken 
at  the  interior  points  z^ ,  ...  z^.  A  similar  divergence  operator  was 

used  by  Peng  [9]  for  effectiveness-factor  calculations  with  polynomial 
basis  functions. 

Application  of  eqns  (20)  and  (21)  to  eqns  (12)  and  (13)  gives  the 
linear  equations 
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HJ  =  L*J  [*iUj)] 


i— Of  •••  n 
j  s  1 ,  • . .  m-1 
k  =  1  f  • • •  n+1 


( z*%\  i  *  ivi  [!tii 

|_  za  dz  dz  jzjj.  J  L  J  Ldz  lzj  _ 

i  =  1 ,  . . .  n 

(23) 

j  —  1  f  •  •  •  n 
k  =  1 ,  •  •  •  n 

from  which  the  matrices  [A  ]  an1  [B  ]  for  the  given  problem  are 

kj  kj 

calculated. 

Insertion  of  eqns  (20)  and  (21 )  into  eqns  (4)-(6)  then  gives  the 
collocation  equations 


- 1  w 


Y_  Vj L  w  • 12  v 


i  =  1 ,  . . .  n  +1 
c 

k  =  1 ,  . . .  n 


Insertion  of  eqn  (24)  into  the  left  side  of  (25)  gives  a  system  of 
equations  for  the  mesh-point  values  ^(z^);  we  solve  this  system  by 
Newton's  method.  The  symmetry  of  the  problems  considered  here  would 
allow  replacing  some  of  eqns  (25)  by  stoichiometric  relations  among 
the  fluxes  [31 ];  however,  the  equations  as  shown  allow  a  simpler 
extension  to  non-symmetric  problems. 
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If  equilibrium  reactions  or  immobile  species  are  included  in  the 
reaction  scheme,  then  a  corresponding  number  of  differential  equations 
must  be  replaced  by  local  equations  as  in  [33 ]•  These  changes  do  not 
arise  in  the  examples  considered  here. 

COMPUTATION  ffiOCEDURE 


A  reference  state  was  chosen  for  each  problem  as  described  in 

Examples  1-3  below.  The  maximum  number  of  independent  production 

rates,  m,  was  computed  as  described  under  eqn  (10).  The  matrices  D(£  ) 

~  ~R 

and.R'C^)  were  computed  analytically  from  the  expressions  for  the 

—1  2 

fluxes  and  reaction  rates,  and  the  eigenvalues  of  [-  D  (£  )  L  R'(C  )] 

~  ~R  ~  ~R 

were  computed  by  the  Q-R  method.  Single  precision  (8  digits)  sufficed 
up  to  this  point. 

The  function  Q  (z)  was  determined  from  eqns  (17)  and  (19)  by  a 
n 

modified  Gram-Schmidt  algorithm.  The  zeroes  of  Q  (z)  were  found  by  a 

n 

grid  search  and  Newton  iteration.  The  weight  matrices  [A^]  and  [E^] 
were  then  ccanputed  from  eqns  (22)  and  (23)  by  UJ  decomposition  [30] 
with  partial  pivoting.  These  calculations  were  done  in  double 
precision  (18  digits).  Polynomial  collocation  was  initiated  in  the 
same  way  except  that  the  eigenvalues  were  replaced  by  zeros. 

Equations  (24)  and  (25)  were  then  set  up  in  single  precision  and 
solved  by  Newton's  method,  starting  from  estimates  ^(z^)  based  on  a 
one-point  solution.  The  examples  that  follow  were  solved  with  a 
general-pirpose  program,  which  does  the  calculations  automatically 
when  the  desired  reaction  model  and  conditions  are  presented  in 
digital  form. 
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EXAMPLE  U  SECOND  ORDER  REACTION  WITH  LARGE  THIELE  MODULUS 

As  our  first  example  we  consider  an  isothermal,  irreversible 
second-order  reaction  in  a  porous  catalytic  slab.  Solutions  sure 
obtained  by  collocation  with  the  hyperbolic  basis  functions,  and 
various  alternative  methods  are  compared. 

The  problem  can  be  stated  in  dimensionless  fora  as  follows, 

=  $2  c*2  -1  <  z  <  1  (26) 

dz2 

c*  =  1  at  z  =  +1  (27) 

and  has  one  independent  rate  of  production;  thus,  m  =  1 .  The  problem 
is  non-linear  and  is  stiff  when  the  Thiele  modulus,  4,  is  large;  it 

provides  a  good  test  of  approximate  solution  methods.  Solutions  are 

*  *.  |  2 

given  here  for  the  dimensionless  boundary  flux,  N  =  dc  /dz  I  ^  =  4  n, 

at  Thiele  moduli  of  10,  100,  and  1000. 

For  these  large  values  of  4,  the  concentration  will  drop  almost 

to  zero  at  the  center  of  the  slab.  Linearizing  the  reaction  term  of 

eqn  (26)  at  the  mean  concentration  of  0.5,  we  obtain  the  single 

2  *2  *  I  2 

eigenvalue  A.  =  -  3(4  c  )/3 c  L  _  =  4  .  The  resulting  basis 
1  |Uo 

functions,  for  j  >  1 ,  are 

♦j  =  (l-z)^-1exp[(z-1  )4]  +  0+z)^_1exp[-(z+1  )4]  (28) 

Several  of  these  functions  are  shown  in  Figure  1 . 

Table  1  shows  the  convergence  of  the  approximations  to  the 
* 

boundary  flux  N  with  increasing  number  of  collocation  points.  The 
hyperbolic  functions  of  eqn  (28)  give  rapid  convergence  throughout. 
With  these  functions,  two-point  collocation  gives  the  boundary  flux 
within  about  1  percent  for  Thiele  moduli  4  >  10. 


TABLE  1 


*  2 

Convergence  of  Solutions  for  Boundary  Flux,  N  =  *  n» 
for  Isothermal  Second-Order  Reaction  in  a  Slab 

!  Number  of  Results  with  hyperbolic  functions  Results  with  splines  [36] 

I  Collocation 

Points  $  =  10  $  =  100  *  =  1000  *  =  1000 


1 

7.9736 

91.356 

969-83 

627.5 

2 

8.0537 

81.127 

814.16 

708.5 

3 

8.1505 

81.689 

817.65 

715.8 

4 

8.1614 

81.662 

816.79 

792.5 

5 

8.1637 

81.661 

816.67 

800.4 

6 

8.1641 

81.654 

816.57 

812.0 

7 

8.1642 

81.652 

816.53 

815.2 

8 

8.1642 

81.651 

816.51 

815.5 

9 

8.1642 

81.650 

816.51 

10 

8.1642 

81.650 

816.50 

11 

8.1642 

81.650 

816.50 

— 

(8.16421 )# 

(81. 6496 )# 

(81 6.497  f 

(816.497)^ 

U 

From  the  general  solution  for  a  single  reaction  in  a  slab  [5,  27,  3]. 


The  exact  concentration  profile  at  $  =  1000  is  shown  in  Figure  2 
along  with  point  values  computed  by  the  present  method.  Good 
agreement  is  found  at  each  n  >  1  except  for  the  point  of  lowest 
concentration,  where  the  linearization  of  the  second-order  kinetics 
is  least  realistic.  Note  that  each  set  of  collocation  points  is 
positioned  nicely  in  the  reaction  zone. 

Villadsen  and  Michelsen  [36]  have  analyzed  this  problem  by- 
collocation  with  global  and  piecewise  polynomials.  Using  global 
polynomials  at  $  =  100,  they  obtained  a  boundary  flux  of  100.0  for 
n  =  8,  and  83.9  for  n  =  12.  From  Table  1  we  see  that  the  hyperbolic 
functions  give  more  accurate  values  than  these,  even  for  n  as  small 
as  2. 

Piecewise  polynomials  (splines)  were  used  by  Villadsen  and 
Michelsen  [36]  to  solve  this  problem  at  $  =  1000,  where  global 
polynomials  failed.  The  general  solution  for  single  reaction  in  a  slab 
[5,  27,  3]  was  used  to  select  the  breakpoints  of  the  piecewise 
polynomials.  Their  best  results  are  shown  in  Table  1;  these  were 
obtained  with  a  breakpoint  at  z  =  0.9967  for  n  >  1 ,  and  another  at 
z  =  0.991  for  n  >  4.  Comparing  the  last  two  columns  of  Table  1,  we  see 
that  each  hyperbolic  solution  (except  for  n  =  1 )  approximates  N  more 
accurately  than  the  2n-point  solution  based  on  splines. 

The  problem  for  $  =  1000  may  also  be  solved  by  a  modified 

Pater son-Cresswell  method.  In  the  original  method  [26],  the  reaction 

is  assumed  to  be  so  fast  that  the  concentration  in  an  interior  region 

[0,z  ]  may  be  neglected,  and  eqn  (5)  need  only  be  integrated  over  the 
s 

region  [z  ,1].  In  the  modified  method  [36]  the  concentration 
gradient,  rather  than  the  concentration  itself,  is  neglected  in  the 
region  [0,z^] ,  and  the  boundary  condition  at  z  =  z^  becomes  dc/dz  =  0 
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Z 


Figure  2.  Concentrations  in  a  catalyst  slab  with  isothermal  second-order 
reaction  at  Thiele  modulus  of  1000.  The  collocation  solutions 
are  based  on  the  functions  of  eqn  (28),  along  with  $q(z)  =  1. 


rather  than  c  =  0.  With  polynomial  collocation,  seme  experimentation 

is  needed  to  select  a  suitable  value  of  z  ;  however,  with  the 

s 

hyperbolic  functions,  the  choice  of  zg  is  not  so  critical. 

For  a  chosen  zg  we  compute  a  solution  with  a  modified  Thiele 

modulus  $  =  (1  -  z  )*.  Let  the  modified  boundary  flux  be  N* ;  then 

ss  s 

.ft 

the  boundary  flux  of  the  original  problem  becomes  N  =  N/(1-z).  T< 

s  s 

obtain  the  boundary  flux  for  *  =  1000,  we  first  try  zg  =  0.9  which 

gives  4>  =  100.  Then  N*  =  81.7  (Table  1,  n  >  3),  and  N*  =  817.;  this 
s  s  - 

is  comparable  to  the  result  obtained  from  our  global  procedure.  The 

alternate  choice  z  =0.99  gives  *  =10,  whence  N*  =  8.16  (Table  1, 
s  s  s 

\  * 

n  >  4),  and  N  =816.  Both  choices  of  z  give  good  approximations  of 
—  s 

* 

the  exact  flux  N  =  816.497;  however,  the  direct  solution  shown  in 
Table  1  is  at  least  as  quick. 

EXAMPLE  2l  THE  WEISZ-HICKS  PROBLEM 

As  a  second  test  of  the  new  collocation  procedure,  we  consider  the 
non- isothermal  effectiveness-factor  problem  of  Weisz  and  Hicks  [37]. 
This  problem  is  described  by  eqns  (4)- (8)  with  spherical  symmetry, 
constant  transport  properties,  and  a  first-order  Arrhenius  kinetic 
model. 

The  kinetic  model  chosen  for  this  example  is  =  -  R^  =  c^ 

exp[-1 5000* (1  /T  -  1/500)].  Here  again,  m  =  1.  The  other  chosen  values 

are  W  =  0,  B  =  0,  D  =  b  ,  L  =  3,  k  =  1 ,  Ap  -  -100,  c  =  1 , 
k  0  ijs  ij  eff  ’0 

and  Tq  =  500  in  eqns  (1 )  and  (3)-(8).  These  values  make  D  a  unit 
matrix,  and  correspond  to  4>  =  3,  S  =  0.2,  and  y  =  30  in  the  notation 
of  Weisz  and  Hicks  [37]. 

For  this  single-reaction  problem,  eqns  (4)-(8)  yield  the  steady 
state  relation  (T  -  Tq)/Tq  =  (1  -  c  )g  throughout  the  particle. 

Equation  (5)  then  reduces  to  the  form  [37] 
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with 


1  L(  z2  dc*  )  =  $2  R*(c*)  (29) 

z^  dz  dz 


# .  *.  # 

R  (c  )  =  c  exp 


Yp(1  -  c*) 

1  +  8(l-e*) 


(30) 


This  function  is  plotted  in  Figure  3  for  the  chosen  values  of  &  and  y 

The  multiplicity  criteria  of  lass  [22],  as  well  as  those  of 

Stewart  and  Villadsen  [32],  predict  a  strongly  ignited  steady  state. 

Thus,  the  reactant  concentration  will  span  the  range  of  Figure  3,  and 

the  reaction  rate  will  have  a  strong  peak  within  the  particle. 

A  linearized  solution  (n  =  m)  will  clearly  not  he  adequate  here; 

the  added  basis  functions  of  eqn  (11 )  will  he  essential.  The  region  to 

the  left  of  the  peak  in  Figure  3  is  chosen  for  the  linearization,  in 

order  to  describe  accurately  the  inner  border  of  the  reaction  zone. 

* 

linearization  at  c  =0  gives  =  1336.,  whereas  an  orthogonal  linear 
fit  of  the  rate  function  in  the  region  of  positive  slope  gives  = 
457.8.  Linearization  at  a  point  to  the  right  of  the  peak  (as  would  be 
appropriate  for  an  unignited  particle)  gives  a  negative,  . lal  A1  and 
yields  polynomial  basis  functions  according  to  the  selection  rules 
described  above. 

Table  2  shows  the  convergence  of  the  solutions  for  the  boundary 

flux.  Both  hyperbolic  function  sets  give  good  results  for  n  >  3;  the 
* 

set  based  on  =  0  is  preferred.  To  obtain  better  than  2  percent 
accuracy  with  this  set  we  need  only  4  collocation  points,  whereas 
with  polynomials  9  points  are  required. 


EXAMPLE  3.  CATALYTIC  REFORMING  OF  C7  HYDROCARBONS 


Krane  and  co-workers  [21 ]  gave  a  kinetic  model  for  catalytic 
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Figure  3.  Steady-state  rate  finction  of  eqn  ( 30)  for  the  Weisz-Hicks 
problem  with  6  =  0.2  and  y  =  30. 
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TABLE  2 

Convergence  of  Solutions  for  Boundary  Flux  for  the 


Weisz-Hicks  Problem  with 

4  =  3.,  6  =  0.2  and  y  =  30. 

Number  of 

Collocation  Solutions, 

*  2 

H  s}  fj/3 

Collocation 

Hyperbolic  functions 

Polynomials 

Points, 

at-  \fT%d1 

at=  \/4577§^ 

v 0 

n 

i  >  0 

i  >  0 

all  i 

1 

33.524 

22.012 

4.942 

2 

12.148 

7.828 

12.955 

3 

8.583 

10.418 

6.661 

4 

9.655 

8.830 

7.317 

5 

9.337 

9.573 

9.218 

6 

9-500 

9-546 

10.359 

7 

9.438 

9-371 

9.130 

8 

9.447 

9-428 

9.209 

9 

9.456 

9-475 

9.546 

10 

9.451 

9.459 

9.502 

11 

9.450 

9.443 

9.406 

12 

9.451 

9.448 

9.445 

13 

9.451 

9.453 

9.465 

a  * 

ft  Determined  by  linearization  at  c  =0. 

■ft.  # 

fflf  Determined  from  a  least-squares  linear  fit  of  R  (c  ) 
over  the  region  of  positive  slope  in  Figure  2. 
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reforming  of  C7  hydrocarbons  over  platinum-alumina.  Their  reaction 
scheme  is  as  follows: 


5 

(0  / 

(s) 

n-Heptane  ^^2^^ 
'\\ 

^  (N) 

4  (A) 

Cracked 

Naphthenes 

c  *  Aromatics 

Products 

"  Isoheptanes 

Their  pseudo-homogeneous  rate  expressions  have  been  generalized  by 
Guertin  et  al.  [13],  with  the  equilibrium  constants  ,  K^,  and 
recalculated  from  the  API  tables  [l]  for  self-consistency .  For  those 
calculations  and  the  following  ones,  species  A  was  taken  as  toluene, 

C  as  an  equimolar  mixture  of  propane  and  i-butane,  N  as 
methylcyclohexane,  and  I  as  an  equimolar  mixture  of  2-  and 
3-methylhexane . 

The  present  theory  requires  local  kinetic  expressions,  rather 

than  pseudo-homogeneous  ones.  To  obtain  these,  we  have  adjusted  the 

rate  constants  to  obtain  mean  rates  <R.>  consistent  with  those  of 

1 

Krane  et  al.  at  the  conditions  of  their  experiments.  A  particle 
diameter  of  1/16"  and  length  of  3/32"  were  assigned  to  those 
experiments  after  consultation  with  John  Sinfelt;  this  gives  an 
effective  sphere  radius  of  0.9  nnn  by  the  Wheeler-Aris  rule  [3,  p.  162]. 
Activation  energies  of  15,000  K  were  assumed  for  the  forward 
reactions;  these  values  are  not  critical  since  the  temperature 
corrections  required  are  small. 

The  transport  parameters  for  a  similar  catalyst  have  been 
determined  by  Feng  [9]  by  fitting  eqn  (1 )  to  his  mass  flux 
experiments.  A  one-point  lumping  of  the  pore  size  distribution  gave 
W1  =  0.081479  and  r^  =  374  A,  whereas  a  two-point  lumping  gave  = 
0.11803,  =  64  A,  W2  =  0.Q32828,  and  =  845  A;  in  both  cases  a 
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c  2  1 

permeability  =  4.5  10  g  cm  sec  atm  was  obtained.  A  thermal 

conductivity  k  =  0.0005  cal  cm  ^  sec  ^  K  ^  is  given  by  Sehr  [29] 
err 

for  a  similar  catalyst;  this  value  was  used  for  the  non-isothermal 

—5 

cases.  A  catalyst  density  of  1 .22  g  cm  J  was  used,  as  for  Peng's 
experiments. 

The  transport  coefficients  y  and  cfi  were  calculated  at  769.24  K 

J 

from  eqns  (1.4-18,19)  and  (16.4-12,15,16)  of  Bird  et  al.  [6];  for 
other  temperatures  a  correction  (T/769.24)  ’  was  applied. 
Lennard-Jones  parameters  were  obtained  for  hydrogen  from  [6],  and 
for  each  hydrocarbon  from  correlation  (iii)  of  Tee  et  al.  [34]  with 
critical  properties  and  vapor  pressures  from  the  API  tables  [l]. 
Surface  diffusion  was  neglected. 

Prom  these  physical  data  and  the  production  rates  <R^>  reported 
by  Krane  et  al.,  we  have  obtained  the  folio*.- rr#  local  :*';e  equations 
for  a  temperature  of  769.24  K  (925  P): 

=  0.702  p-°*63  [  ps  -  Pj/2.95  ] 

R,  =  0.567  p'1  [Pg-PH  p/1 .77  ] 

=  0.567  p"1  [  pT  -  PH  Pn(2.95/1  .77)  ] 

(31) 

£  =  78.8  p"1  [  pw  -  p3  p  /377700.  ] 

4  N  ri  A 

=  0.106  p”0,67  pg 

Jt  *  0.061  p-0,67  p 

6  I 

With  these  local  expressions,  the  computed  mean  production  rates  <R^> 

•  **+  >»■ . - 
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closely  approximate  the  values  reported  in  [21 ]  for  the  n-heptane  and 
methylcyclohexane  conversion  experiments. 

The  pseudo-homogeneous  rate  constants  of  Krane  et  al.  [21 ,  13] 
differ  from  our  local  values  by  the  following  factors:  0.52  for 
reaction  1;  0.63  for  reactions  2  and  3;  0.13  for  reaction  4;  1 .47  for 
reactions  5  and  6.  These  comparisons  show  that  the  intraparticle 
gradients  are  important,  even  for  this  small  particle  size.  Similar 
conclusions  were  reached  by  Hettinger  et  al.  [14]  from  their 
experiments  with  several  particle  sizes. 

Calculations  are  reported  here  for  catalytic  reforming  in 
spherical  particles  with  the  following  boundary  conditions  at  z  =  1 : 


Pg  =  PN  =  PA  =  1  atm,  pH  =  15  atm,  p];  =  pc  =  0,  T  =  769-24  K.  The 
large  naphthene  partial  pressure  corresponds  to  a  first-reactor  inlet 
condition;  the  fate  of  the  naphthenes  in  that  region  is  crucial  to  the 
selectivity  of  the  process.  Calculations  are  given  for  spherical 
particles  with  radii  of  0.9  and  1 .6  mm;  the  latter  size  corresponds 
approximately  to  a  1/8”  x  1/8"  pelleted  catalyst. 

The  reference  state  was  chosen  here  as  ?(z^ )  from  a  one-point 
polynomial  collocation  solution;  other  choices  of  would  give 
comparable  results  for  this  mildly  non-linear  reaction  model.  Analysis 
of  the  stoichiometry  by  the  method  of  [31  ]  gives  m  =  4  as  the  maximum 
number  of  independent  production  rates. 

Table  3  shows  the  values  and  for  the  0.9  mm  particle  with  a 

one-point  pore  size  distribution  model.  The  complex  eigenvalues  arise 
from  the  temperature  dependence  of  the  reaction  rates;  the  eigenvalues 
were  all  real  in  the  isothermal  case. 


Table  4  shows  the  convergence  of  the  boundary  fluxes  for  two 
species  with  increasing  n;  the  efficiency  of  the  hyperbolic  functions 
is  evident,  especially  for  hydrogen.  The  fluxes  of  both  species  are 
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TABLE  3 


Selection  of  Basis  Functions  for  Catalytic  Reforming 

in  a  Nonisothermal  Particle.  L  =  0.9  nn  and  n  =  1 

w 

Eigenvalues 

Real(VT7) 

k 

Selection 

k 

X 

s  a 

Order(f) 

k 

k,  initial 

after  aQ  =  < 

(*) 

-0.2  10-7 

(•) 

(•) 

0.3  10"7 

(») 

(») 

-0.00119  +  0.00443  i 

(»> 

1 

-0.00119  -  0.00443  i 

0.0474 

4,  8,  ... 

2 

17.381  +  0.56871  i 

4.1696 

3»  7t  • • • 

3 

17.381  -  0.56871  i 

4.1696 

2,  6,  ... 

4 

422.57 

20.556 

1 »  5,  ... 

!*)  Since  this  problem  has  m  s  4, 

the  4  largest  values  Re(VoT) 

are  selected.  The  other  eigenvalues  are  not  used. 

(ff)  These  are  not  the  indices  j  of  eqn  (14);  the  latter  are 
obtained  by  labeling  the  chosen  a  values  in  ascending  order 

as  a  ...  a  .  The  grouping  operations  are  then  performed. 

0  n 


* 


TABLE  4 


Convergence  Test  of  Boundary  Fluxes  for  Catalytic  Reforming 

in  a  Nonisothermal  Particle.  L  =  0.9  ma  and  n  =  1 

w 


Number  of 
Collocation 
Points 

L  N  ,  mole  cm' 
iz 

-1  sec"1  * 

106 

Hyperbolic  functions 

Hydrogen  n-Heptane 
flux  flux 

Polynomial  functions 

Hydrogen  n-Heptane 
flux  flux 

1 

1.4274 

-.06692 

.3709 

-.05133 

2 

1.4279 

-.03992 

.8670 

-.04834 

3 

1.4343 

-.04035 

1.2421 

-.04102 

4 

1.4350 

-.04122 

1.3945 

-.03984 

5 

1.4309 

-.04299 

1.4293 

-.04105 

6 

1.4297 

-.04322 

1.4322 

-.04797 

7 

1.4286 

-.04337 

1.4304 

-.04294 

8 

1.4279 

-.04347 

1.4290 

-.04329 

9 

1.4279 

-.04351 

1.4282 

-.04344 
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accurate  within  1#  for  n  >  6  with  the  hyperbolic  functions,  versus 
n  >  8  with  polynomials.  Since  the  number  of  operations  required 
to  solve  eqn  (26)  varies  essentially  as  n^,  the  hyperbolic  functions 
give  practical  accuracy  in  about  half  the  computing  time. 

Table  5  shows  the  mean  production  rates  computed  with  n  =  7,  for 
several  conditions  and  several  versions  of  the  transport  equations. 

The  first  two  lines  give  results  for  a  nonisothermal  0.9  mm  particle 

with  n  =  1  and  2  in  eqn  (1).  The  appreciable  difference  between  these 

w 

two  cases  is  caused  by  the  wide  distribution  of  pore  sizes  in  the 

catalyst.  The  next  two  lines  are  obtained  for  an  isothermal  particle 

with  n  =  1  and  2.  As  expected,  there  are  only  moderate  differences 
w 

between  the  non-isothermal  and  isothermal  solutions  because  the  net 

heat  of  reaction  is  small.  Line  5  represents  a  simplified  model  in 

which  all  off-diagonal  elements  in  the  F(r)  matrix  are  neglected,  thus 

giving  each  molar  flux  in  the  form  N,  =  -  D  dc./dz.  The  resulting 

iz  ieff  i 

approximations  for  <R^>  are  good  within  10  percent;  larger  deviations 
may  be  expected  for  systems  in  which  the  gaseous  phase  is  not  so 
diluted  with  a  highly  mobile  component. 

Line  6  gives  rates  calculated  with  a  permeability  twice  the 
measured  value.  These  results  agree  closely  with  line  2;  thus, 
approximate  values  of  and  y  should  normally  suffice. 

Lines  7  and  8  give  the  mean  production  rates  for  larger  and 
smaller  particles:  1 .6  mm  and  0  mm.  The  rates,  and  their  ratios, 
vary  strongly  with  particle  size;  the  rate  for  n-heptane  actually 
changes  sign.  Thus,  a  pseudo-homogeneous  kinetic  model  is  not 
appropriate,  at  least  for  these  reaction  conditions. 

Tables  6  and  7  show  the  concentration  and  temperature  profiles 

for  both  particle  sizes,  calculated  from  the  full  model  with  n  =2. 

w 


TABLE  5 

Mean  Production  Rates  for  Catalytic  Reforming  in  Spherical  Particles, 

Computed  with  n  =  7 


I 

i 


<\>, 

mole  g”1 

sec 

,05 

Specifications 

L,  mm 

n-Heptane 

Iso¬ 
heptanes  Naphthenes  Hydrogen 

Toluene 

Cracked 

Products 

Non iso thermal , 

n  =1 
w 

0.9 

-13.17 

44.59 

-191.5 

433.7 

156.9 

6.91 

Non iso thermal, 

n  =2 

w 

0.9 

-13.44 

47.04 

-202.7 

458.2 

165.8 

6.98 

Isothermal , 

n  =1 
w 

0.9 

-13.96 

44.71 

-193.5 

440.6 

159.2 

7.16 

Isothermal , 

n  =2 
w 

0.9 

-14.36 

47.23 

-205.0 

466.2 

168.6 

7.26 

Non iso thermal, 

D  modelC*) , 

ieff 

n  =2 
w 

0.9 

-13.81 

48.80 

-213.1 

480.9 

173.1 

7.20 

Nonisothermal, 

B  doubled, 

0 

n  =2 
w 

0.9 

-13.40 

47.03 

-202.6 

453.6 

165.8 

6.98 

Nonisothermal, 

s3 

II 

l\> 

1.6 

-9.43 

27.53 

-115.5 

264.6 

95.0 

5.80 

Small-particle  limit 

0.0 

29.59 

250.3 

-1489. 

3327. 

1205. 

8.49 

(*)  The  D  model  is  obtained  by  neglecting  the  right-hand  term  of  eqn  (2b). 
ieff 


1 
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TABLE  6 

Profiles  for  Catalytic  Reforming  in  a  Spherical  Particle  of  Radius  0.9  mm 

Computed  with  n  =  7  and  =  2 


!  - - — - - —  . . . 

j  ,  g 

!  Concentrations,  mole  cm-5  *  10° 


Iso—  Cracked  Temp, 

z  n-Heptane  heptanes  Naphthenes  Hydrogen  Toluene  Products  K 


TABLE  7 

Profiles  for  Catalytic  Reforming  in  a  Spherical  Particle  of  Radius  1 .6  am 

Computed  with  n  *  7  and  =  2 


7  6 

Concentrations,  mole  cm-^  *  10 


z 

n-Heptane 

Iso¬ 

heptanes 

Naphthenes 

Hydrogen 

Toluene 

Cracked 

Products 

Temp. 

K 

1.0000 

15.84 

0. 

15.84 

237.6 

15.84 

0. 

769-24 

0.9913 

15.25 

1.02 

11.70 

238.5 

18.80 

0.16 

768.51 

0.9540 

13.00 

3-74 

3.41 

240.3 

24.83 

0.80 

767.01 

0.8854 

9-88 

6.00 

0.63 

240.9 

27.08 

1.84 

766.44 

0.7851 

6.93 

7.48 

0.32 

240.9 

27.60 

3.13 

766.33 

0.6438 

4.80 

8.13 

0.31 

240.8 

27.84 

4.49 

766.29 

0.4641 

3.55 

8.15 

0.30 

240.7 

27.97 

5.73 

766.28 

0.2451 

2.96 

7.92 

0.30 

240.6 

28.00 

6.63 

766.29 
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The  naphthene  and  n-heptane  concentrations  decrease  strongly  toward 
the  center  of  the  particle,  whereas  the  toluene,  isoheptane  and 
cracked  product  concentrations  increase.  The  total  intraparticle 
temperature  drop  is  about  3  degrees  Kelvin  for  both  particle  sizes. 

CONCLUSION 

A  generalized  model  and  numerical  procedure  have  been  presented 
here  for  steady-state  simulation  of  porous  catalyst  particles.  The 
examples  show  the  ability  of  the  method  to  deal  with  fast  or  slow 
kinetics,  multicomponent  diffusion  and  multiple  reactions. 

The  new  collocation  method,  with  basis  functions  given  by  eqns 
(14),  allows  efficient  solutions  both  for  stiff  and  non-stiff 
differential  equations.  For  non-stiff  equations  (small  eigenvalues) 
the  method  reduces  to  a  global  polynomial  collocation  scheme.  For 
stiff  equations  a  similar  reduction  eventually  occurs  at  large  values 
of  n,  through  the  grouping  rules  described  above  eqns  (14). 

The  selection  of  the  Oj values  still  requires  some  judgment. 
Research  on  this  aspect  of  the  method  is  continuing. 

Example  3,  as  well  as  the  data  in  [14],  show  that  intraparticle 
gradients  are  important  in  catalytic  reforming  operations. 

Use  of  two  pore  sizes  in  eqn  (1 )  is  recommended  for  wide  pore 
size  distributions,  such  as  commonly  occur  in  reforming  catalysts. 
This  is  more  accurate  than  the  assunption  of  a  single  pore  size,  and 
takes  very  little  more  computing  time. 
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NOTATION 


weights  for  dimensionless  gradient  operator,  eqn  (20) 

geometric  parameter,  0  for  slabs  and  2  for  spheres 

adjustable  coefficients,  eqn  (12),  consistent  units 

-2  -1 

permeability  coefficient,  g  cm  sec  atm 
adjustable  coefficients,  eqn  (IT) 
column  vector  with  elements  c^ 

_3 

total  molar  concentration  in  pore  space,  mole  can 
molar  concentration  of  species  i  in  pore  space,  mole  cm  ^ 


JB  (e,T) 

V', 


~s 

E 

kj 

F(r) 

Ah 

keff 


M. 


=  c  /c  ,  dimensionless  concentration  of  reactant  species 

1  I  f 

matrix  with  elements  D 
row  i  of  D 

transport  coefficient  calculated  from  eqns  (A.2)-(A.5)» 
consistent  units 

2  -1 

binary  diffusivity  of  pair  ij,  cm  sec 

=  r\/  52  R  T/9  ,  Khudsen  diffusivity  of  gas  i  in  a  pore 

2  -1 

of  radius  r,  can  sec 

2  -1 

surface  diffusion  coefficient  matrix,  cm  sec 

adjustable  coefficients,  eqn  (15),  consistent  units 

weights  for  dimensionless  divergence  operator,  eqn  (21 ) 

—2 

matrix  with  elements  defined  by  eqn  (2),  sec  cm 
partial  enthalpy  of  species  i,  cal  mole-1 
enthalpy  of  reaction,  cal  mole-1 

effective  thermal  conductivity  of  porous  medium,  cal  cm-1  sec-1  K' 
half-thickness  or  radius  of  particle,  cm 
mass  number  of  species  1,  g  mole-1 
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m 


N 

iz 


maximum  number  of  independent  production  rates  Ri  permitted  by 

stoichiometry  and  local  constraints;  see  discussion  under  eqn  (10) 

column  vector  with  elements  N  ,  . . .  N 

~1  ~nc+1 

column  vector  with  elements  N  ,  . . .  N 


~1 


~n„ 


-2  -1 

smoothed  flux  of  species  i,  mole  cm  sec  ,  for  i  =  1 

—2  —1 

smoothed  total  energy  flux,  cal  cm  sec  ,  for  i  =  E 


n 


z-component  of 


R’ 

R 


r 

T 


2 

=  4>  n/(a*1  )*  dimensionless  reactant  flux  at  z  =  1 

number  of  interior  collocation  points 

number  of  gaseous  species 

number  of  pore  sizes  in  eqn  (1 )  - 

total  pressure,  atm 

column  vector  with  elements  pi 

partial  pressure  of  species  i  in  pore  space,  atm 

leading  factor  in  residual  functions,  eqn  (1 6) 

specific  rate  of  reaction  j,  mole  g-1  sec-1 

column  vector  with  elements  R^ 

=  p  )  v..R.  ,  rate  of  production  of  species  i, 

P  L-.  J1  J 

i 

mole  cm-^  sec”^ ,  for  i  =  1 ,  . . .  n  ;  R.  =  0  for  i  =  E 

C  1 

average  of  R^  over  the  particle  volume 

matrix  with  elements  3Ri/3Cj,  consistent  units 

gas  constant,  consistent  units 

pore  radius,  A 

temperature,  K 

porosity-tortuosity  coefficient  in  eqn  (1 ),  dimensionless 
*  c^/c,  mole  fraction  of  species  i  in  pore  space 


z  dimensionless  coordinate  relative  to  L,  measured  from  center 

of  particle 

Greek  letters 


8 

Y 

6 

id 

n 

»k 

u 

Vji 

5 

5i 

pp 

* 

♦4 


dimensionless  parameter  in  basis  function  eqn  (14) 

dimensionless  heat  of  reaction  in  Ref.  [37] 

dimensionless  activation  energy  in  Ref.  [37] 

Kronecker  symbol,  unity  when  i  =  d  and  zero  otherwise 

=  <R  >/R  n,  effectiveness  factor  for  single  reaction 

dimensionless  eigenvalues  of  eqn  (9) 

viscosity  of  gas  mixture,  g  cm-1  sec-1 

stoichiometric  coefficient  of  species  i  in  reaction  j 

column  vector  with  elements  ^ 

=  c.  for  i  =  1 ,  ...  n  ;  S.  =  T  for  i  =  E 
i  c  i 

particle  density,  gam 

Thiele  modulus 

basis  function,  eqn  (14) 


Subscripts 

E  energy;  element  n  +1  in  column  vectors  N  and  £. 

c  ~ 

H  molecular  hydrogen 

R  reference  state  for  1 inearization 

s  auxiliary  variables  in  Pater son-Cresswell  method 

z  z  component 

0  outer  surface 

Superscripts 


per  unit  mass  of  catalyst 

approximation  in  terms  of  basis  functions;  see  eqns  (12)  and  (13) 
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APPENDIX  A.  THE  MATRIX  D 

T 

Insertion  of  the  gas-law  expressions  g  =  cRT  and  p  =  1  g  into 
eqn  (1 )  (here  is  a  row  vector  with  all  elements  unity)  gives 
nw 

N(m)  =  -V  W  [P(r jr^Jc  +  c^lnT]  -  D  Wc 
k=1 

-  (BqRT/ y)  [c  ij®]  Jc  -  (B0cR/y)  c  (A.1) 

This  gives  the  following  elements  of  D: 


For  i  /  E  and  j  ^  E, 


3Ni 


aVc 


MM  tel 


l\ tP(rk)r;  ♦  (B0W/„)  ct  ♦  »  (*.2) 


For  i  f  E, 


die  = 


3N.  A 


3Vt 


=  V  V.  [F(r  )-1c]  /T  +  (B  dR/y)c.  (A.3) 

?  L k  • k  ‘ 1  0  1 


The  remaining  elements  of  D  are  obtained  from  eqns  (3)  and  (A. 2): 
For  d  +  E, 


i=1 


H, 


(A. 4) 


APPENDIX  B.  PROPERTIES  OP  THE  GRID  POINTS  AND  BASIS  FUNCTIONS 


Let  a  .  ...  a  be  distinct  real  numbers,  and  let  P_(z),  ...  P  (z) 
Os  0  s 

be  polynomials  of  degrees  m^,  ...  mg.  Then  the  function  ^P^(z)  exp(a^ 


z) 


has  no  more  than  -1  + 


Y_  ^mi+  1 ) 


zeros  on  the  real  axis,  unless  it 


*  r  n 

vanishes  identically;  see  Polya  and  Szego  [28J,  p.  48,  Theorem  75*  This 

result  can  be  used  to  bound  the  number  of  zeros  of  any  function 

M 

Fm(z)  =  ^  cj  ♦jt*)  (B.1) 

j=0 

constructed  from  the  first  M  +  1  members  of  the  set  in  eqns  (14). 

For  the  slab,  eqn  (B.1 )  can  be  written  in  the  form 

PM(z)  =  P0(1_Z)  +P0(l4z) 


Pi(1-z)  exp[-o^(1-z)] 


iX> 


Pj^O+z)  exp[-c^(1+z)] 


(B.2) 


iX) 


by  use  of  eqns  (14).  Here  the  are  the  distinct  members  of  the  set 
{“o’  •••  sum  +  1)  is  (n^  +  1)  +  2(1^  +  1)  +  ...  ■ 


2M  +  1 ;  hence,  by  the  above  theorem  of  Polya  and  Szego,  the  number  of 

zeros  of  F  (z)  on  the  real  axis  does  not  exceed  2M.  Since  F  (z)  is  evei 
M  M 

the  number  of  zeros  on  the  real  interval  (0,<»)  then  does  not  exceed  M. 


j 
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For  the  apnere,  eqns  (B.1)  and  (11)  give 
z  Fm(z)  =  z  PQ(1-z)  -  z  PQ(1  +z) 


y  P^l-z)  exp[-ot^(l-z)] 


iX) 


-  y  Pi(l+z)  expC-a^l+z)] 


(B.3) 


i>0 


The  sum 


Y  <*. 
u  1 


+  1 )  is  now  2M  +  2;  this  allows,  at  most,  2M  +  1  real 


zeros  for  z  F„(z),  and  one  less  for  F  (z).  Again,  F  (z)  is  even  and 
M  MW 

has  no  more  than  M  zeros  on  the  real  interval  (0,°°). 

Now  let  v  he  the  number  of  sign  changes  of  Qr(z)  on  the  real 
interval  (0,1 )  in  a  slab  or  sphere.  The  result  just  shown  ensures  that 
v  <  n;  it  also  ensures  that  there  exists  an  expansion  F^(z)  given  by 
(B.1 )  whose  sign  agrees  on  (0,1 )  with  that  of  ^(z) .  Than  the  integral 
/1q  (z)  F  (z)  z9^1 dz  is  nonzero,  since  its  integrand  is  positive 

0  n  v 

except  at  the  zeros  of  Q^(z) .  This  result  is  compatible  with  eqns  (19) 

and  (B.1 )  if  and  only  if  v  >  n.  But  v  <  n  ,  as  noted  above.  Hence, 

v  =  n;  that  is,  Q  (z)  has  exactly  n  sign  changes  on  (0,1 ).  Bach  sign 
n 

change  is  a  zero,  since  Q  (z)  is  continuous. 

n 


The  determinant 


M+1 


VV 

4  (z  ) 
1  1 


W 


•  •  • 

•  •  • 

•  •  • 


4>  (z  ) 

1  M+1 


Vvi} 


(B.4) 


is  nonzero  as  long  as  the  real  numbers  z^ ,  ...  z^  are  distinct  and 


positive,  since  no  linear  combination  of  $rt(z),  ...  4>„(z)  can  vanish 

U  M 

at  more  than  M  distinct  points  on  (0,»).  Therefore,  the  matrix 


[«J>  (z  )]  of  eqn  (22)  (for  which  M  =  n)  is  non-singular,  and  the  weight 

*“  J 

matrix  [A^]  is  uniquely  defined.  The  uniqueness  of  [E^]  in  eqn  (23) 
can  be  proved  similarly,  using  the  basis  set  {d<j>..(z)/dz} . 
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